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Abstract 

We present a method of general apphcabihty for finding exact or accurate approximations to 
bond percolation thresholds for a wide class of lattices. To every lattice we sytematically associate 
a polynomial, the root of which in [0, 1] is the conjectured critical point. The method makes the 
correct prediction for every exactly solved problem, and comparison with numerical results shows 
that it is very close, but not exact, for many others. We focus primarily on the Archimedean 
lattices, in which all vertices are equivalent, but this restriction is not crucial. Some results we find 
are kagome: Pc = 0.524430..., (3, 12^) : = 0.740423..., (3^ 4^) : p^ = 0.419615..., (3, 4, 6, 4) : p^ = 
0.524821..., (4,8^) : = 0.676835..., (3^,4,3,4): = 0.414120... . The results are generally within 
10^^ of numerical estimates. For the inhomogeneous checkerboard and bowtie lattices, errors in 
the formulas (if they are not exact) are less than 10~^. 
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I. INTRODUCTION 

Since its introduction in the 1950s by Broadbent and Hammersley [l, 2|, percolation has 
been a rich source of interesting and challenging problems. Its very simple construction has 
made it the classic model of disordered media and phase transitions. However, despite its 
apparent simplicity, exact answers to even the most basic questions have been difficult to 
come by. 

The bond percolation process is defined on a lattice by declaring each edge to be open 
with probability p and closed with probability 1 —p, where p G [0, 1]. One can then ask many 
questions about the resulting connected bonds, such as the average size of open clusters, the 
cluster density, and the probability of crossing large regions. We are concerned here with 
the critical threshold, pc, defined as the probability at which an infinite connected cluster 
first appears. Percolation thresholds depend on the lattice under consideration, with each 
one presenting a different set of challenges. However, there is a limited number of cases in 
which the critical probability is known exactly, with all currently solved problems confined 
to a particular self-dual class of lattices, as discussed below. Although we will be concerned 
mainly with the bond problem here, we could also consider site percolation in which p is 
the occupation probability of the vertices (sites) on the graph. There are likewise very few 
exactly known site thresholds. 

In the absence of exact solutions, we must seek alternative approaches. One is to study the 
problem numerically, and thresholds, both site and bond, are now known to high precision for 



a wide variety of systems 



have also been found 
and Walters 



am 
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3. 



Mathematically rigorous results in the form of inequalities 



17H191, and the bounds are continually being narrowed. Riordan 



20| have developed rigorous "99.9999% confidence intervals," subtly distinct 



from rigorous bounds, in which simulations are used to place percolation thresholds in a 
certain range with very high probability. Another avenue of pursuit that has emerged is the 
search for "universal" formulas; closed expressions for the percolation threshold as a function 
of a few lattice parameters, usually focussing on coordination number and dimension (see 



2l| and [22| for just two examples). These formulas generally meet with mixed success, 
providing good approximations for some lattices but performing poorly for others. Whether 
or not one believes a simple closed formula is a reasonable goal, it is clear that the next 
best thing is to have a well-defined procedure for finding a good analytical approximation to 



2 



the threshold of any given lattice. We present such a procedure here for bond percolation. 
With each lattice we show how to associate a polynomial, the solution of which in [0, 1] is 
the prediction for the critical point. Some of the main results of this work were reported 
in a recent Letter 23 1. Our approximation is related to Wu's homogeneity assumption 



for the Potts model as both methods give similar results when applied to the same lattices. 
Wu's work goes back several decades, but recently he has applied it to a variety of new 



lattice systems 



25|, though not the ones presented here. 



Although our method works for any two-dimensional graph with a periodic structure, 
we confine our attention here to the Archimedean lattices, which are composed of regular 
polygons such that all vertices are equivalent. There are 11 such graphs and they are shown 
in Figure [H where we have used the Griinbaum and Shephard |26j notation: given any 
vertex, one lists, in order, the number of sides of the polygons adjacent to it. For example, 
the kagome lattice (Figure [T]i), is (3,6,3,6). If the same number would appear more than 
once in succession, an exponent is used, so the lattice in Figure [1^ is named (3^,4^) rather 
than (3,3,3,4,4). 



II. METHOD 



Currently, the most general tool available for deriving exact percolation thresholds is the 



triangle-triangle transformation 



27 



, which is based upon duality and is a generalization 



of the star-triangle transformation 281]. The method works only for lattices in which the 
basic cell is contained within three sites and arranged in a self-dual way 29|. The simplest 
example of a self-dual array of triangles is just the triangular array itself, as shown in Figure 
[21 Solutions to these problems are derived by considering probabilities of events that take 
place on a single cell of the given lattice. For a lattice of^the type in Figure [2l the critical 
threshold is identified using the condition given by 
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3Q|: 



PiA,B,C) = PiA,B,C) 



(1) 



where P{A, B,C) is the probability that the three vertices can be connected through open 
bonds in the cell, P{A, B, C) is the probability that none are connected, and the shaded 



a) 



b) 



c) 



d) 




i) 



j) 



k) 



FIG. 1: The Archimedean lattices; a) square; b) honeycomb; c) triangular; d) kagome; e) (4,8^); 
f) (3, 122); g) (33^ 42). (3, 4, 6, 4); i) (3^, 4, 3, 4); j) (3^, 6); k) (4, 6, 12) . 





a) 



b) 



FIG. 2: a) A class of exactly solvable lattices; b) every shaded triangle of a) can represent any 
network of sites and bonds contained in the vertices {A, B,C). 



triangle can represent any network of sites and bonds. This condition is also related to Wu's 
criticality condition qA = C (where A and C are elements of the Boltzmann factor of a 
single triangle) for the g— state Potts model on triangular arrays 31|. 

If the cell is a simple star (Figure |3]) where each bond is given different probabilities, 



4 



A 



B 



FIG. 3: Using this simple star in Figure [2)3 results in the honeycomb lattice (Figure Hb)- 

p,r,s, then ([T]) gives the criticality condition 

H{p, r, s) = prs — rp — rs — ps + 1 = 



(2) 



which is the inhomogeneous threshold for the honeycomb lattice. By assigning a different 
probability to each bond of the unit cell, we have obtained a critical surface. We can find 
the homogeneous critical probability, pc, by setting all probabilities equal, giving 

H{p,p,p) =p^ -3p^ + 1 = (3) 



so pc= 1 -2sin7r/18 ^ 0.652704 [28|. 



A useful property of two-dimensional bond thresholds is the duality relationship. For a 
given graph, L, its dual, Ld, is constructed as shown in Figure HI a vertex is placed in every 
face of L and connected by bonds to all the vertices in neighbouring faces. Bond thresholds 
of a dual pair of lattices are related by 28 1 



Pc{L) = l-p,{Ld). (4) 

This can be generalized, for lattices obeying ([!]), to inhomogeneous thresholds. The dual of 
the honeycomb lattice is the triangular lattice, so we can write 

T(p, r, s) = — H(l — 1 — r, 1 — s) (5) 

giving 

T(p, r, s) = prs — p — r — s + l= (6) 

which is the result we would derive using ([T]) directly. The homogeneous polynomial is 
p^ ~ 3p + 1 = with solution pc = 2sin'7r/18 ~ 0.347296 consistent with (jlj). The minus 
sign in ([5]) is arbitrary given that we set these functions equal to 0, but it ensures that the 
constant term in the critical function is +1, a convention we will always employ. 



FIG. 4: The duality transformation. The honeycomb and triangular lattices are dual pairs. 

It is clear that any inhomogeneous problem, with probabilities {pi,P2, ■■■,Pn) on the bonds 
of the unit cell, that can be solved with ([T]) has a threshold of the form 



where no probability appears as a power greater than 1, a property we will refer to as 
"linearity" . The central assumption of the method presented here is that a// bond percolation 
thresholds have this form. The goal is then to find the function / for each lattice. This will 
be done by demanding that the function reduce correctly to all known special cases. For 
example, in the honeycomb lattice case, setting p = 1 contracts the p bond to zero length 
resulting in a square lattice. Equation ([21), predicts S(r, s) = 11(1, r,s) = 1— r — s = 0, 
which is the correct threshold |28|. The challenge is to find sufficient conditions for each / 
to completely constrain the function. 

One might worry that for lattices that are not handled by ([1]), the requirements of linearity 
of / in all its probabilities and its correct reduction in all special cases might sometimes be in 
conflict. Interestingly, this does not appear ever to be the case, i.e., such an / always seems 
to exist. Another potential difficulty is that there may be more than one / that satisfies all 
the conditions and thus the method will give more than one answer. However, we will show 
that if / exists, it is unique. 




(7) 
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a) 



b) 



FIG. 5: a) The martini lattice; b) The assignment of probabihties on the unit cell. 




a) 




FIG. 6: a) The martini-A lattice; b) The assignment of probabilities on the unit cell. 

III. EXACT SOLUTIONS 

In what follows, we will need several exact solutions in order to have cases to which our 
approximations must reduce. We list a few of the important ones here. 





a) 



b) 



FIG. 7: a) The martini-B lattice; b) The assignment of probabilities on the unit cell. 



A. Martini lattice 

We start with the martini lattice 
inhomogeneous bond threshold 
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33| (Figure [5^), and recover the previously reported 
34l |. To solve this problem, we could use ([T]) with Figure 
\5jp, which requires a rather involved enumeration of paths on the unit cell. However, by 
virtue of ([1]), and the fact that the probabilites P are all linear functions of the bond 
probabilities, the critical function, denoted M{pi,p2,P3,ri,r2,r3), is also linear in all its 
arguments. Furthermore, we know that setting ri = results in a honeycomb lattice in 
which two of the bonds are doubled in series, i.e., 

M(pi,p2,P3,0,r2,r3) = H(pi,p2r3,P3'"2)- 

In order to satisfy this condition and linearity, the function must be of the form 

^{Pi,P2,P3,ri,r2,r3) = H(pi,p2'^3,P3'^2) +ri(l){puP2,P3,r2,r3) (8) 

where is as-yet undetermined. Setting ri = 1 gives another honeycomb lattice, this time 
with a more complicated vertical bond (Fig [8]). The probability of crossing this bond is 
Pi[l — (1 — ?"2)(1 — ''"3)], i-e., the pi bond must be open and the r2 and bonds cannot both 
be closed. Thus, 

M(pi,P2,P3,l,'^2,'^3) = H(pi[l - (1 -r2)(l -r3)],p2,P3)- (9) 
Using ([9]) in ([8]) allows us to determine (j) uniquely. The resulting expression for M is 
^iPi,P2,P3,ri,r2,r3) = (1 - ri)H(j9i, J92r3,p3r2) +riH(]9i[l - (1 -r2)(l -r3)],p2,P3) 



and expanding this out leads to the final result [3l|, 



^iPl,P2,P3,ri,r2,r3) = 1 - PlP2'^3 - ^2^3^ - PlP3'^2 - PlP2n^2 

- Pmnrs - p2P3r2r3 + PiP2Pzrir2 

+ PmPsnrs + piP2P3r2rs + piP2rir2r3 

+ PiP3rir2r3 + p2P3rir2r3 - 2piP2P3'"i'"2'"3 = . (10) 

Solving the polynomial M{p,p,p,p,p,p) = gives the exact homogeneous bond percolation 
threshold, Pc = 1/ V^- Two special cases will also be useful later; the martini-A lattice 
(Figure |6]) found by contracting the ps bond, MA(pi,p2 5 ''"i, "'^2, ''^3) = M(pi,p25 1, '"i, "'^2, ^"3), 



Pi 




FIG. 8: The unit cell resulting from setting ri = 1 in the martini lattice. This is just a honeycomb 
lattice with a complicated vertical bond. 

and the martini-B lattice (Figure [7]) which results from contraction of the p2 and bonds, 
MB(pi, ri, r2, rg) = M(pi, 1, 1, ri, ra, rg). 

The method we used to find this solution is somewhat less work than using ([T]), however 
its greatest virtue lies in the fact that it can be generalized to lattices that are not handled 
by ([T]). As we will show, the thresholds predicted in those cases are not exact but are always 
very good approximations, usually holding to 4 or more decimal places. 



B. Stack of three triangles 

The stack of three triangles (Figure |9]) will be useful in deriving the approximation for 
the (3,4,6,4) lattice. The unit cell has 9 probabilities, and we denote its critical function 
by SoT(pi,p25P35 ''^i, "'^2, '"s, Si, S2, S3) . This is a very long formula, and it is included as a 
text file among the source files of this submission. Here, we just report the polynomial for 
the homogeneous case: 

SoT(p, p, p, p, p, p, p,p,p) = l- - 9p^ + 3p^ + 45p^ - 72p^ + 45p'^ - 12p^ + p^ = . 



The critical point of this lattice is Pc = 0.471629... 35|. 



C. Rocket lattice 



This is another lattice that is helpful for the (3,4,6,4) approximation. It is shown in 
Figure [TOl which we call "rocket" from the shape of the unit cell. The critical function is 
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a) b) 
FIG. 9: a) The stack of three triangles; b) the assignment of probabihties. 




a) b) 

FIG. 10: a) The "rocket" lattice; b) the assignment of probabilities. 

denoted R(p, r, s, t, m, f , w, x), and once again we will forego reporting the full threshold but 
have included it in the supplemental material. The result for the homogeneous case is, 

1 - - Ap^ + 13/ - I2p^ + 3/ = . 

The critical probability is Pc = 0.669182... . 

D. Decorated square 

This lattice (Figure [TTl) will be useful when we consider the (3^, 4, 3, 4) lattice. It is really 
just a square lattice with a complicated horizontal bond. Its threshold is straightforward to 
find: 

DS(p, r,s,t,u,v) = 1 — p — st — ru — rtv — suv + rstv 

+ rstu + rsuv + rtuv + stuv — 2rstuv = 0. (11) 
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a) 



b) 



FIG. 11: a) The "decorated" square lattice; b) the assignment of probabihties. 

s 




r 

I 



FIG. 12: The assignment of probabihties on the kagome lattice unit cell. 



IV. APPROXIMATIONS TO UNSOLVED PROBLEMS 

We turn now to lattices whose thresholds are not known exactly. To derive these results 
we use the same reasoning as above, but because none of these cases is handled by ([T]), 
we would not expect any of them to be exact. This expectation is generally confirmed by 
comparison with precise numerical results. However, as we will show, the approximations 
they provide are consistently very good, and two interesting cases may turn out to be exact. 



A. Kagome lattice 



Bond percolation on the kagome lattice (Figure [ffi) has been the subject of particular 



curiosity and numerous conjectures over the years [2^, 129|, |36|]. This is the one unknown 
threshold for the four most basic lattices that are generally studied (the square, triangular, 
honeycomb, and kagome). The attraction of this problem also stems from its similarity to 



the double-bond honeycomb lattice, which has an exactly known bond threshold 



28|. Most 
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of the attempted solutions have been based on finding a simple connection between the two, 
but none have been successful (although some have been quite close). 

A shortcut like the one we used for the martini lattice is possible here. We assign the 
probabilities as shown in Figure [T21 If we set the si bond probability to the result should 
be the martini-A lattice, so we must have 



K(ri,r2,r3,Si,S2,S3) = MA(s3, S2, r2, ra, ri) + Si0(ri, r2, rs, S2, S3). (12) 

Setting this bond to 1 gives the martini-B lattice with its pi bond doubled in parallel. We 
should therefore set pi = 1 — (1 — S2)(l — S3) in the B lattice threshold: 

K(ri,r2,r3,l,S2,S3) = MB(1 - (1 - S2)(l - S3), n, ra, r3). (13) 

Setting si = 1 in (fT2|) and substituting ( !T3|) allows determination of 0(ri, r2, r3, S2, S3), giving 

K(ri, r2, r3, Si, S2, S3) = (1 - Si)MA(s3, S2, r2, rs, ri) + SiMB(l - (1 - S2)(l - S3), ri, r2, r^). 

(14) 

After expanding, the result is: 

K(ri, r2, r3, Si, S2, S3) = 1 - riSi - r2S2 - r3S3 - Sir2r3 - S2rir3 - S3rir2 

- ?^iS2S3 - r2SiS3 - r3SiS2 + Sirir2r3 + S2rir2r3 
+ S3rir2r3 + rir2SiS3 + rir3SiS2 + r2r3SiS2 

+ '"2'"3-5lS3 + rir2S2S3 + r2SiS2S3 + r3SiS2S3 

+ rir3S2S3 + riSiS2S3 - rir2r3SiS3 - rir2r3S2S3 

- rir2r3SiS2 - rir2SiS2S3 - rir3SiS2S3 

- ?^2?^3SlS2S3 + rir2r3SiS2S3 . (15) 



This inhomogeneous threshold was reported previously in [34'], but was derived by a different 
method. The critical probability is the solution to K{p,p,p,p,p,p) = 0, i.e., 

1 - 3p2 - 6p^ + 12/ - + / = (16) 

with solution in [0, 1], Pc = 0.524430..., which is identical to the value conjectured by Wu in 
:.979 |2J]. Two recent numerical estimates arep^ = 0.52440499(2) [Sj andpc = 0.52440510(5) 
1^, so we are close, but not quite right. The inhomogeneous threshold (JTSil is remarkable 
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FIG. 13: Assignment of probabilities on the (4,8^) lattice. 



in that, absent a high-precision numerical result (which has only been available for the last 
twelve years |37:]), it would be easy to convince oneself that it is the correct answer. It has 
all the necessary symmetry properties, and makes no incorrect prediction for any special 
case one can consider. Moreover, as is obvious from the derivation, it is the only function 
linear in its arguments that could possibly serve as the kagome threshold. The fact that it 
is incorrect means that the linearity assumption is wrong, but it is interesting that we can 
find such a function regardless, and that it makes a prediction so close to the correct answer. 
These properties are characteristic of all the criticality functions we will find here. 

Equation ( |T4|) has the formal appearance of an average of two lattice thresholds. This 
is the general form of a linear critical surface. If we have a lattice L\ with critical surface 
fiiipi, ■■■,Pn) = 0, that reduces to L2, which is a lattice with m < n bonds in its unit cell 
and threshold /l2(^1) •■•^'"■m), when pi = 1 and L3, which has k < n bonds in its unit cell 
and threshold fisisi, Sk), when pi = 0, then we have 



where pi and rji are functions that are at most first-order in their arguments. Li is a kind 
of "average" of L2 and L3. Equation f|T71) is the only way one can form a threshold that 
is linear in probabilities and that correctly reduces to L2 and L3 in the appropriate limits. 
This indicates that, in general, if a linear function exists that makes correct predictions in all 
special cases, then it is the only such function and the method of obtaining it is irrelevant. 



/li(Pi, = PlfL2{Pl{P2, ■■■,Pn), Pm{P2, ■■■,Pn)) 




(17) 
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B. (4,8^) lattice 



The basic cell of the (4, 8^) ("four-eight" or "FE" lattice) lattice has six bonds and is 
shown in Figure US] along with the probability assignments. Removal of the t bond results 
in the honeycomb lattice with two doubled bonds, and setting t = 1 gives the martini-A 
lattice. Equation (fT7|) implies 

FE(p, r, s, t, u, f ) = (1 — t)H(p, ru, sv) + tMA{u, v, r, s,p). 

Expanding gives the complete function, which will also be useful later, 

FE(p, r, s, t,u,v) = 1 — {pru + stu + psv + rtv) — {rsuv + ptuv) + 
prsuv + ptruv + pstuv + rstuv + prstu + prstv — 2prstuv = 0. 

The critical polynomial FE{p,p,p,p,p,p) = is then 

1 - 4p^ - 2p^ + 6p^ -2p^ = 

with solution on [0,1], Pc = 0.676835... . According to Parviainen 4|], this lattice has 
Pc = 0.676802... with standard error 6.3 x 10"^. Although our answer lies outside the error 
bars, it does agree to four significant figures. It also lies well within the confidence interval 
of Riordan and Walters 20|] shown in Table [B 

This threshold makes another interesting prediction. Setting u = v = 1 gives the square 
lattice with four different probabilities as shown in Figure [HI and we have 

CB(p, r,s,t) = 1 — pr — ps — rs ~ pt — rt ~ st + prs + prt + pst + rst = 0. (18) 



This arrangement was named the "checkerboard" by Wu [38| and he conjectured its threshold 
for the g— state Potts model based on symmetry and special cases, in much the same way as 
we are doing here. Our result is identical with the o — )■ 1 limit of Wu's formula. Although 



his conjecture was shown to be incorrect for q = 3 [39|, our simulations suggest that it may 
be exact for percolation. However, proving this will require a new technique, as it is not 
handled by 1^. 



14 



p p 





V t 

p 




s 

t 

p 




s 

r 

P 




s 

r t 




s 




s 



FIG. 14: The square lattice with the checkerboard assignment of probabihties. 




FIG. 15: The assignment of probabilities on the (3, 12^) unit cell. 
C. (3, 122) lattice 



The approximation for the (3, 12^) lattice (often called the 3-12 lattice 40|]) was also 
previously mentioned by the authors [3^. We will just repeat the result here, which we 
write in a factorized form to save space, 

1 - mim2(r3 + rir2-rir2r3)(s3 + SiS2-SiS2S3) 

- mim3(r2 + rir3 - rir2r3)(s2 + S1S3 - S1S2S3) 

- m2m3(ri + r2r3 - rir2r3)(si + S2S3 - S1S2S3) 
+ mim2ms{rir2 + rir^ + r2r3 - 2rir2r3) 

X (siS2 + S1S3 + S2S3 - 2S1S2S3) = , (19) 
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a) b) 

FIG. 16: a) The "striped" triangular lattice; b) "re-partioned" version of a). The critical threshold 
should be invariant under this permutation of the probabilities. 

where we have used the probability assignments from Figure [151 The polynomial for the 
homogeneous case factorizes to 

(1 + p - + p^)(l - p + p2 + - + Ap') = . 

The part that contains the solution in [0, 1] is 

+ + -7p^ + 4p^ = (20) 



and yields Pc = 0.7404233179... . According to Ding et al. 
so our approximation agrees to 5 significant figures. 



16|, Pc(3,122) ^ 0.7404207(2), 



D. (3^42) lattice 

For the remaining lattices, we must employ a more brute-force method to find the critical 
function, requiring the use of Mathematica to perform the algebra. For the (3^,4^) lattice, 
we consider an inhomogeneous triangular lattice, which we call the "striped" triangular 
lattice because the probabilities are assigned in distinct stripes (Figure fT6k). The (3^,4^) 
lattice appears when ri = 0. We start by defining the most general function that is at most 
first-order in its variables, 

1 1 

f{pi,ri,si,p2,r2,S2) = ^...^aijMmnPlris'lp^r^s^ (21) 

i=0 n=0 
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where there are 64 a's to be determined (in general, there will be 2" coefficients, where n is the 
number of bonds in a unit cell). However, we are free to set the constant term to 1, because 
we can multiply / = by any constant and not change the threshold. Now we impose all 
the known symmetries and special cases onto /. For example, simply re-partitioning the 
lattice as shown in Figure [TBb should not change the result, so 

/(pi,ri,Si,p2,r2,S2) = /(pi,ri,S2,P2,r2,Si). 

This constraint tells us that a certain set of the a's are equal to each other, reducing the 
number of coefficients by 16 in this case. Similarly, the threshold should be invariant under 
a flip about a vertical axis: 

/(Pl,n,Sl,P2,r2,S2) = /(ri,Pi,Si,r2,P2,S2) 

which eliminates another 18 coefficients. Next, there is the symmetry between the 1 and 2 
bonds, 

/(pi,ri,Si,p2,r2,S2) = /(P2,r2,s2,pi,ri,si) 

eliminating 9 a's, leaving 20 to be determined. Setting ri = si = should give the threshold 
for the martini-B lattice, 

/(pi,0,0,p2,r2,S2) = MB(pi,S2,P2,r2) 

leaving 9 undetermined coefficients. Because we are dealing with the triangular lattice 
Pi = P2 = P,fi = r2 = r, Si = S2 = s should yield (j6]). The function / will not simply 
reduce to the function T(p, r, s) because f{p,r,s,p,r,s) has second-order terms. Rather, 
T(j9, r, s) must factor out of the resulting expression. Because T(p, r, s) = is equivalent to 

p = (r + s — 1)/ {rs — 1) we must have 

/r + s - 1 r + s -1 \ 

f r, s, — , r,s] =0. 

\ rs — 1 rs — 1 J 

This leaves only one coefficient unconstrained. We find this by observing that ri = r2 = 
gives the same lattice as setting ri = p2 = 0, with some probabilities changed, i.e., 

/(pi,0,Si,P2,0,S2) = /(pi,0,Si,0,P2,S2). 
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This finally fixes the threshold, which is given by 



ST(pi, ri, si, p2, ^^2, S2) = 1 - si - S2 - P1P2 - rir2 - pin - - ^2^2 - Pir2 

+ S1S2 + PiP2r2 + PiP2ri + Pirir2 + P2rir2 + PiViSi 

+ P2r2Sl + PiriS2 + P2r2S2 - PiriSiS2 - P2r2SlS2 

- PiP2rir2S2 - PiP2rir2Si + PiP2rir2SiS2 = (22) 
The (3'^, 4^) lattice is obtained by setting ri = 0: 

TF{pi,p2,r2,Si,S2) = 1 - Si - S2 + S1S2 - P1P2 - Pir2 - P2r2 

+ P2r2Si + P2r2S2 + PiP2r2 - P2r2SiS2 = (23) 

The homogeneous polynomial is 

TF(p, p, p,p,p) = 1 — 2p — 2p^ + 2)p^ — = 

with solution on [0,1] pc = 0.419308... . Parviainen finds pc ~ 0.41964191 with standard 
error 4.3 x 10^'', a difference of 0.00033. Once again, our result is close but not in agreement 
with the numerical value. It may seem somewhat surprising that we could satisfy all our 
various criteria using only a function first-order in probabilities. What is more, we have a 
formula that makes predictions for cases that we did not explicitly consider in the derivation. 



For example, setting si = gives an inhomogeneous version of the bow-tie lattice 41| . 
BT(j),r, s,t,u) = ST(t, s,0,r,p,u) (Figure fTTI) . so 

BT(p, r,s,t,u) = 1 — u — pr — ps — rs — pt 
— rt — st + prs + prt + pst 

+ rst + pru + stu — prstu = 0. (24) 

Setting the probabilities equal gives 

l-p-6p'^ + 6p^ -p^ = (25) 



with pbowtie ^ 0.404518..., which is the exact answer found by Wierman [4l|]. Based upon 
numerical studies, we believe that the more general equation (!24|) may also be exact 42| . 
even though it is not obvious how to derive it using the normal duality arguments. 
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Another interesting prediction is found by setting ri = r2 = (Fig [18]). This again results 
in a square lattice, but the probabilities are assigned differently from the checkerboard. We 
call this the "striped" square, and the threshold is predicted by (12^ to be 

- P1P2 + S1S2 - si - S2 + 1 = 0. (26) 

This formula has several interesting properties. For example, if we set Pi = P2 = P and 
■Si = S2 = s we have the standard square lattice and fl26|) can be factored to give 

(1 +J9- s)(l -p- s) = 

leading to the correct threshold 1 — p — s = S{p, s) = 0. We can discard the first term in 
brackets because s = 1+p obviously does not give answers in (0, 1). Although the underlying 
lattice is self-dual, this assignment of probabilities is not. The dual lattice is again a square, 
but the inequivalent p and s bonds are interchanged. Our formula reflects this property too, 
because if we set Pi ^ 1 — Pi and Sj — )■ 1 — we indeed end up with 

- S1S2 + P1P2 - Pi - P2 + 1 = 0. 

Finally, setting all probabilities equal leads to = 1/2 as required for the square lattice. 
Although this formula rather nicely exhibits all the required features, it does not turn out 
to be exact. This is easily seen by a numerical example. If we set Si = S2 = P2 = 0.4, then 
predicts pi = 0.9. However, our numerical work j42| gives a value of pi ~ 0.9013, which 
rules out this prediction. 

Finally, we point out that the path we took to constrain ST(pi, ri, Si,p2, r2, S2) might not 
be the most efficient or direct; the order in which we apply the conditions is mostly just the 
order in which we thought of them. However, as we pointed out above, once one has found 
a consistent linear threshold it is the only function that has all the correct properties. The 
converse is not true however — two different lattices may be self-consistently represented 
by the same linear function. An example is equation f l2B]) . which is also the exact threshold 
of the square lattice with vertical bonds replaced by two in series, and horizontal bonds 
replaced by two in parallel (Figure (fT9|) ). i.e. S{piP2, 1 — (1 — si)(l — S2)). 
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FIG. 17: The 5-bond bow-tie lattice. The predicted threshold for this situation, equation (i24|) . is 
possibly exact. 
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FIG. 18: The lattice resulting from setting ri = r2 = in (j22p with the unit cell highlighted. 
Equation ([26]) is the predicted threshold. 




FIG. 19: Unit cell of the square lattice for which ()26p is the exact threshold. 



20 




FIG. 20: Probabilities for the (3, 4, 6, 4) lattice. 



E. (3,4,6,4) lattice 

The basic cell for the (3,4,6,4) lattice of Figure [T]i is shown in Figure [201 It has 12 
bonds, making it the largest cell we will consider. There are 4095 coefficients that must 
be constrained (the constant is set to 1 as usual). The constraints that fully determine the 
critical function are as follows: 

1. We can rotate the unit cell by 120 degrees without changing the threshold, meaning 
that we have 

f{puP2, P3, P4, ri,r2, rg, r^, Si, Sa, S3, S4) = /(rg, rg, ri^r^, 53, S3, Si, S4, p2, Ps, Pi,Pa). 
Rotating the opposite direction does not give any additional constraint. 

2. It is possible to re-partition the lattice, as we did in Figure [T6] for the striped triangular 
situation, so the unit cell is still a set of three squares but the probabilities are shuffled, 
i.e., 

/(pi,P2,P3,P4,n,?^2,r'3,'^4,Sl,S2,S3,S4) = / (Sl , S3, S4, S2, ^4, P2, Pi , PS, ^^2 , ^4, r3, Ti) . 

3. We can reflect the lattice about a vertical line drawn through the middle of the unit 
cell, so 

/(Pi,P2,P3,P4, ri, r2, rs, r4, si, S2, S3, S4) = /(p2,Pi,P3,P4, S2, Si, S3, S4, r2, ri, r3, ta). 
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Reflections about other axes do not give any further constraints. 

4. Contracting the p^, r^, S4 bonds, i.e., setting their probabihties to 1, resuhs in the stack 
of triangles, so 

f{pi,P2,P3, 1, n, ^2, rg, 1, Si, S2, S3, 1) = SoT(pi, ps, PS, ^2, Tg, Si, S2, S3). 

5. Setting S3,r3 = results in the rocket lattice (Figure [TOl) with a comphcated v bond: 
/(Pi,P2,P3,P4,ri,r2,0,r4, Si,S2,0, S4) = R(r4, S4,p4,pi,p2, 1- (1-P3)(l-S2ri), Si, r2)) 

6. If we set r2 = = S2 = S3 = 0, the (4, 8^) lattice with some doubled bonds results. 
Because we previously found the only linear function of probabilities that would qualify 
as that threshold, (|T8l) . the current function must reduce to that one. That is, 

/(Pi,P2,P3,P4, n, 0, rs, 0, si, 0, 0, S4) = FE(p4,pi,p2,P3, S1S4, rirs) 

7. A version of the kagome lattice results when we remove the ps, ri, and S2 bonds. Once 
again, the only consistent linear threshold is given by f|T5|) so we expect reduction to 
that formula: 

f{Pi,P2, 0,P4, 0, r2, rs, r4, si, 0, S3, S4) = K(p4, S4, r4, S3r3,p2r2, sipi). 

8. Setting r2 = rs = Si = S3 = 1 gives a (3^,4^) with several decorated bonds. This 
requires 

fiPl,P2,P3,P4,ri, 1, l,r4, 1, S2, 1, S4) = 
TF(1 - (1 - p,){l - p2),P3, 1 - (1 - ri)(l - r4), 1 - (1 - S2)(l - S4),P4). 

9. Setting Pi = P2 = ^3 = Si = 1 and S4 = gives a version of the dual to the (4, 8^) 
lattice. As such we will need 

FED(p, r, s, t, M, v) = -FE(1 - p, 1 - r, 1 - s, 1 - t, 1 - m, 1 - i;). 

and so 

/(I, 1, P3, P4, n, r2, 1, r4, 1, S2, Ss, S4) = 

FED(S4,1 - (1 -P3)(l -P4),l - (1 -r2)(l - Ss),ri,S2,r4) 
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10. We can get a version of the honeycomb lattice by assigning pi = ri = r2 = S3 = S4 = 0: 

/(O, P2, P3, P4, 0, 0, rs, r4, Si, S2, 0, 0) = H(p2P3P4, S2r-3) 



11. We now have only one coefficient left. Setting ps = Pi = results in a 10-bond 
version of the dual to the bow-tie lattice, with the known threshold Pc = 0.595482... . 
Requiring 

/(Pc, Pc, 0, 0, Pc, Pc, Pc, Pc, Pc, Pc, Pc, Pc) = (27) 

fixes the final coefficient. 

All coefficients are now constrained, and the full expression, containing 879 terms, is included 
in the supplemental material. The polynomial for this lattice is 

1 - - 12/ - 6p^ + 69p^ + 60p^ - 363p® + 448^^ - 252p^° + 66p^^ - 6p^^ = 

with Pc = 0.524821111889... . Parviainen's numerical result is pc = 0.52483258 with standard 
error 5.3 x 10"^ Once again we are outside the error bars but the difference is only 0.000012 

As with the (3^,4^) lattice, we have a few other predictions of this formula that are 
worth discussing. As mentioned in step 9, setting ps = p^ = gives a 10-probability 
version of the bow-tie dual lattice, from which the 10-probability bow-tie (Fig [22|) threshold, 
EBT{pi,ri, Si,ti,ui,p2,r2, S2,t2,U2), can be found. As usual, all special cases seem to have 
the right properties. To test whether this inhomogeneous threshold is actually exact, we will 
reduce it to the extended checkerboard, i.e., the square lattice with 8 different probabilities 
(Fig[25h). by setting pi = p2 = I and p^ = Pi = 0. So 

ECB(pi, ri, si, ti,p2, r2, S2, ^2) = /(1, 1, 0, 0,p2, ^2, ^2, t2,pi, si, ti, n). 

Now, setting ri = and the other probabilities equal, we find a lattice with an unknown 
threshold (Fig [23b). The polynomial is 

1 - 3p^ - 5p^ + I2p^ - 7p^ + = 

with solution on [0, 1] Pc = 0.5696764123... . Numerical simulation of this lattice gives Pc ~ 



0.56982 



42], so unfortunately the extended bow-tie and checkerboard thresholds predicted 



here cannot be exact. 
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FIG. 22: The bow-tie lattice with 10 different probabilities. 
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FIG. 23: a) The extended checkerboard lattice with 8 independent probabilities; b) The result of 
setting ri = 0. We use this lattice to test our prediction for a). 
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FIG. 24: a) The line graph of the square lattice, which is used to derive the (3^, 4, 3, 4) polynomial; 
b) The assignment of probabilities. 



F. (3^,4,3,4) lattice 

The final case we consider is the (3^,4,3,4) lattice. To do this, we start with the lattice 
shown in Figure which is the covering graph of the square lattice with 12 probabilities. 
We will constrain it until we have a sufficient number of coefficients for the lattice we want, 
which appears when si = t2 = 0. An approximation that is helpful here is the one for the 
4-bond square matching (SM) lattice (Figure [25]) . By now it should be clear how this is 
done, so we only state the result: 

+ P1P2S2 + P1S1S2 + P2S1S2 - 2piP2SiS2 = (28) 
For uniform probabilities we have 

SM(p, p,p,p) = l-4p + Ap^ - 2/ = 
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with solution pc = 0.266385..., which of course is not exact |42 |. 

We will now simply list the constraints used in deriving the approximation to the 
(32,4,3,4) lattice. 

1. Rotation: 

/(Pl,P2,P3,P4,n,?^2,?^3,'^4,Sl,S2,tl,t2) = / (P2 , P4, Pi , Ps, ^^2, ?^4, ^ , ^^3, S2 , Si , ^2, ^l) 

2. Reflection: 

/(Pl,P2,P3,P4,?"l,?"2,?^3,?^4, Sl,S2,tl,t2) = / (?^2, ^1 , Tg, , Pi , P4, P3, ^1 , ^2, , ^2) 

3. Re-partitioning: 

/(Pi, P2, P3, P4, ri, r2, rs, r4, Si, S2, ti, t2) = /(n, r2, rg, r4,Pl,P2,P3,P4, tl, ^2, Si, S2) 



f{Pi,P2,P3, P4, ri, r2, rs, r4, Si, 0, h, 0) 
EBT(pi,p2,P3,P4,ri,r2,r3,r4,Si,ti) 



5. 



/(1,P2,P3,P4, I,r2,r3,r4, Si,S2,tl,t2) = 

ST(p^, r^, r4, s*, t*, P4) 

where = 1 - (1 - ti)(l - rg), = 1 - (1 - r2)(l - t2), = 1 - (1 - Si)(l - pg), 
and = 1 - {1 - P2){1 - S2) . 

6. Using ([II]): 

/(pi,P2,P3,0,ri,r2,r3, 1,0,0,^1,^2) = 
DS(pi,p2,P3,ri, 1 - (1 - ti)(l - r2), 1 - (1 - t2)(l - ra)) 



7. Using (|23]): 



f{Pl,P2,P3,P4,ri, 1, l,r4, Si,0,ti,t2) = 

TF(1 - (1 - ri)(l - ti)(l - t2)(l - r4),P3, 1 - (1 - Pi)(l - P4), ^1,^2) 
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8. Now we need ([28D: 

/(Pl,P2,P3,P4, 1, 1, 1,^4, Si,S2,tl,t2) = 

SM(1 - (1 - p2)(l - Ps), 1 - (1 - Pi)(l - Pa), Su S2)) 

9. 

/(Pl,P2,P3,P4,ri,r2, l,r4, Si,S2, 1, 1) = 

SM(1 - (1 - p2){l - Pa), 1 - (1 - Pi)(l - P4), Si, S2) 

10. 

/(Pl,P2,P3,P4, l,r2,r3,l,Si,S2,tl,l) = 
SM(1 - (1 - P2)(l - Ps), 1 - (1 - Pl)(l - P4), Si, S2) 

11. 

/(Pl,P2,P3,P4,l,'^2, 1,^^4,51,52,^1,1) = 

SM(1 - (1 - p2)(l - Pa), 1 - (1 - Pi)(l - P4), Si, s,) 

12. 

/(Pl,P2,P3,P4,n,?^2,?^3,'^4,0,l,l,0) = CB(p,f,S,t) 

where p = 1 — (1 — r2)(l — r4), f = 1 — (1 — P3)(l — P4), s = 1 — (1 — Pi)(l — P2), and 
t = l-(l-ri)(l-r3). 

This finally constrains the coefficients we need to get the critical surface for the (3^,4,3,4) 
lattice, which has 458 separate terms and can be found in the supplemental material. The 
homogeneous polynomial is 

1 - 4p2 - 12p3 - 2p^ + 106p^ - 186p^ + 132p^ - 36p*^ - 2p^ + 2p^° = 

or pc = 0.414120304... . Parviainen's numerical result is Pc ~ 0.41413743, with standard 
error 4.6 x 10"^. Although outside the error bars, our prediction differs by only 0.000017. 
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a) 




b) 



FIG. 25: a) The square matching graph. The crossed bonds are not connected inside the square; 
b) The assignment of probabihties. 



lattice 


approx 
Pc 


„num 
Pc 


bounds 


|^approx_^num 


kagome 


0.5244297 


0.5244050 


[0.52415,0.52465] 


2.5 X 10"^ 


(3,122) 


0.7404233 


0.7404207 


[0.7402, 0.7407] 


2.6 X 10-6 


(4,82) 


0.676835 


0.676802 


[0.6766,0.6770] 


3.3 X 10-5 


(3^42) 


0.419615 


0.419642 


[0.4194,0.4199] 


2.7 X 10-5 


(3,4,6,4) 


0.524821 


0.524833 


[0.5246,0.5251] 


1.2 X 10-5 


(32,4,3,4) 


0.414120 


0.414137 


[0.4139,0.4144] 


1.7 X 10-5 



TABLE I: Comparison of bond percolation estimates with numerical results (the kagome result is 



from Feng, Deng and Blote 



Parviainen 



a, 



(3, 12 ) is from Ding, Fu, Guo and Wu [I6l ]. and the rest are from 



and the Riordan and Walters confidence intervals 



0- 



G. Extensions 

One question that now arises is whether there is any benefit to considering larger unit cells 
for any of our predictions. The extended checkerboard (Figure [23]) provides one opportunity 
to study this issue. For example, we can see what prediction the extended checkerboard 
formula makes for the simple checkerboard by setting Pi = P2 = Pi ri = r2 = r, Si = S2 = 
s, ti = t2 = t. The resulting threshold can be written in factored form: 

[1 — pr — ps — rs + prs — pt — rt + prt — st + pst + rst) x 
{1+pr+ps — rs— prs — pt + rt — prt + st — pst — rst + 2prst) = 

We recognize the first term in brackets as our original prediction for the checkerboard, (fTSll . 
and the second term can be discarded as it does not have a root in [0, 1]. Evidently, we get the 
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same prediction for the checkerboard, even if we use a larger unit celL However, if we consider 
again the striped square lattice shown in Figure [HI the story is different. Our extended 
checkerboard formula makes a prediction for this case as we can get Figure [TS] from Figure 
[23] by substituting pi S2,ri pi,Si Pi,P2 Si,r2 -> P2,S2 ^2,^2 P2- 

This prediction is 

1 - 2pip2 + pIpI - si- 2p\siS2 - SP1P2S1S2 
+ 8pIp2SiS2 - 2plsiS2 + 8piplsiS2 

- QpIpIsiS2 + 2plsls2 + 'ipiP2sls2 

- ^pIp2sIs2 + 2plsls2 - 4piP2S?S2 

+ 2p\p\s\s2 - sl + 2p\sisl + AP1P2S1S2 

- 4:pIp2SisI + 2plsisl - Apiplsisl 
+ '^pIpIsisI + slsl - 2pls\sl 

- 2pip2slsl + 2plp2s\sl - 2plsjsl 

+ 2piplslsl = (29) 

and equation ( I26l) does not factor out. Now, if we set si = S2 = P2 = 0.4, we find pi = 
0.901308, which, compared to the previous answer pi = 0.9, is much closer to the numerical 



value of pc = 0.901263(2) 42|]. Clearly some refinement of the estimate is possible by 
considering larger unit cells, at least in some cases. 

Similarly, we can improve the estimate for the (3^,4^) lattice, by considering the trian- 
gular lattice with 12 different probabilities (Figure [26]) rather than 6 as for the striped 
triangular problem (Figure [T6i) . We will again list only the constraints used to find 

f{pi, ri, Si,P2, r2, S2,P3,?^3, S3,P4,r4, S4). 

1. Re-partitioning: 

/(pi,ri,Si,P2,'^2,S2,P3,'^3, 53,^4,^^4,54) = /(P2,'^2,S2,Pl,ri,Si,P4,r4,S4,P3,r3,S3) 

2. Re-partitioning: 

f{pi,ri,Si,p2,r2,S2,P3,r3,S3,p4,ri,Si) = /(p3, ^^3, ^3,^4, ^"4, S4,Pi, n, Si,p2, ^2, S2) 

3. Re-partitioning: 

/(pi,ri,Si,p2,r2, 52,^3,^^3, 53,^4,^4,54) = /(P2,r4, 53,pi,r3,S4,P4,r2,Si,P3,ri,S2) 
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4. Reflection: 

/(Pl,n,Si,P2,r-2,S2,P3,'^3,S3,P4,r-4, S4) = /(S4, r4, ^4, S2, P2, S3, Ts, ps, Si, Ti, pi) 

5. Reflection: 

/(pi,ri,si,p2,r2, 52,^3,^^3, 53,^4,^4,54) = /(s2,ri,P3,S4,r3,pi,Si,r2,P4, S3,r4,p2) 

6. Using the 10-bond bow-tie from Figure (l22l) : 

f{pi, ri, Si,P2, 0, S2,P3, 0, S3,P4, r4, S4) = EBT(pi, Si,P3, S2, ri,P4, S4,P2, S3, ^4) 

7. 

/(pi, ri, si, p2, ^"2, S2, 0, rg, S3, 0, r4, S4) = EBT(s4, rg, Si, ri, pi, S3, r4, S2, r2, P2) 
8. Using ([22]): 

/(pi, ri, si, p2, r2, S2, 0, 0, 0, p4, r4, S4) = ST(s2, P2, r2, Si, pi, ri, S4, P4, ^-4) 

9. 

/(pi,ri,Si,p2,r2,S2,P3,r3, l,P4,r4, 1) = 
ST(1 - (1 - pi)(l - r,){l - ps), ri, Si, 1 - (1 - p4)(l - r4)(l - P2), r^, S2) 

10. Finally: 

/(pi, 1, Si,P2,r2, S2,P3, 1, S3,P4,r-4, S4) = 
ST(1 - (1 -pi)(l - S4)(l - Si),P4,r4, 1 - (1 - S3)(l -P3)(l - S2),P2,r2) . 

The homogeneous (3^,4^) threshold is given by f{p,0,p,p,p,p,p,0,p,p,p,p) = 0: 

1 - 4p2 - 12p^ + 104p^ - 193p^ + 146p^ - 45p^ + 2p^° = (30) 

for which the solution on [0, 1] gives Pc = 0.419615..., the value given in Table [H This is an 
improvement on the 6-bond estimate of Pc = 0.419308... compared to Parviainen's numerical 
estimate of Pe = 0.4196419(4). 

Finally, we mention that we also found the 12-bond kagome threshold. However, this 
makes the same prediction as the 6-bond case, which we know is incorrect. Evidently, to 
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FIG. 26: Probability assignments on the extended striped triangle. 



get any improvement on this threshold would require a large number of probabilities. We 
have seen a few examples where thresholds are unchanged by employing larger unit cells. 
All those results are either exact or, in the case of the checkerboard and 5-bond bow-tie, 
thought to be exact. However, we can see here that this property cannot be taken as a test 
for exactness, as the kagome threshold exhibits some cell-size independence as well. 



V. DISCUSSION 

The results are summarized in Table [H along with the numerical results of Parviainen 
4|, Feng, Deng and Blote jsl, and Ding, Fu, Guo, and Wu 16 1 as well as the confidence 
intervals of Riordan and Walters {20]. The final two Archimedean graphs, the (3^,6) and 
(4, 6, 12) lattices have large unit cells, containing 15 and 18 bonds respectively. This presents 
a challenge to the method and we have not yet found the polynomials for those systems. 
However, for the Archimedean solutions we did find, the largest difference between our 
predicted value and the numerical result is 0.000033, for the (4, 8^) lattice, and none are 
ruled out by the Riordan and Walters confidence intervals. 

What is not clear, however, is why this linearity argument works as well as it does. If 
we were to give every bond of a lattice a different probability, then the critical function of 
this situation would certainly be "linear" in this infinite number of probabilities. This is 
because the threshold could be located (in principle anyway) by enumerating open paths 
that span increasingly large regions. In any of these paths, a given bond would only be 
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mentioned once, and therefore only the first power of its probabihty would appear in the 
solution. Evidently, by extending the inhomogeneous probabilities over ever greater regions, 
we can approach the exact threshold, and our approximation in which the probabilities cover 
only a single unit cell of the lattice is the first step in this procedure. However, because 
the solutions we found are so close to the numerical values, we have effectively shown that 
convergence to the exact threshold is very fast. 

Naturally, we would like to learn how to go beyond these lowest-order results. Although 
we saw examples of this when we found fl29|) as a refinement of (!26|) and improved the (3'^, 4^) 
threshold, finding the higher-order approximations to most problems will likely prove very 
difficult using the method we have presented here. In our derivations we have been able to 
avoid the question of what our critical functions actually represent, i.e., perhaps they arise 
from the application of some comparison of probabilities like condition ([1]), so it is conceivable 
that there is a more computationally straightforward way to find them. Discovering such a 
specification would be valuable in extending the approximation. 
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